library(foreign)
library(sandwich)
#library(tikzDevice)
library(Matching)
require(causalsens)
data(lalonde)
psid <- read.dta("psid_controls.dta")
psid$u74 <- 1 * (psid$re74 == 0)
psid$u75 <- 1 * (psid$re75 == 0)
psid <- psid[,c(3:ncol(psid), 2)]
names(lalonde) <- names(psid)

mean(lalonde$re78[lalonde$treat == 0]) - mean(psid$re78)

nonexp <- rbind(subset(lalonde, treat == 1), psid)

sens.form <- as.formula(paste("y.adj ~ treat +age + education +
  black + hispanic + married +
  nodegree + re74 + re75 + u74 + u75"))
eff.form <- as.formula(paste("re78 ~ treat+age + education +
  black + hispanic + married +
  nodegree + re74 + re75 + u74 + u75"))
ps.form <- as.formula(paste("treat ~ age + education +
  black + hispanic + married +
  nodegree + re74 + re75 + u74 + u75"))


pmodel <- glm(ps.form, data = lalonde, family = binomial())
ymodel <- lm(eff.form, data = lalonde)
alpha <- seq(-4500, 4500, by = 250)

ll.sens <- causalsens(ymodel, pmodel, ~ age + education, data = lalonde,
                      alpha = alpha, confound = one.sided.att)
plot(ll.sens, type = "raw")
plot(ll.sens, type = "r.squared")


##tikz(file = "../figs/lalonde.tex", width = 5.5, height = 3)
par(mfrow = c(1,2), cex.main = .8, cex.axis = 0.8, cex.lab = 0.8,
    mar = c(4.2,4,1,1) + .1)
plot(x = ll.sens$sens$alpha, y = ll.sens$sens$estimate,
     ylim = c(min(ll.sens$sens$lower), max(ll.sens$sens$upper)), type = "l",
     xlab = "Amount of confounding (α)", ylab = "Estimated effect", bty = "n",
     yaxt = "n", xaxt = "n")
axis(side = 1)
axis(side = 2, las = 2)
abline(h = 0, col = "grey")
abline(v = 0, col = "grey")
abline(h = ll.sens$sens$estimate[ll.sens$sens$alpha == 0] - 1000, lty = 2, col = "black")

polygon(x = c(ll.sens$sens$alpha, rev(ll.sens$sens$alpha)),
        y = c(ll.sens$sens$lower, rev(ll.sens$sens$upper)),
        col = rgb(0.5, 0.5, 0.5, alpha = 0.5),
        border = NA)

plot(x = sign(ll.sens$sens$alpha)*ll.sens$sens$rsqs,
     y = ll.sens$sens$estimate,
     ylim = c(min(ll.sens$sens$lower), max(ll.sens$sens$upper)), type = "l",
     xlab = "Variance explained by confounding \n($R_α^2$)",
     ylab = "Estimated effect", bty = "n", yaxt = "n", xaxt = "n")
axis(side = 1)
axis(side = 2, las = 2)
abline(h = 0, col = "grey")
abline(v = 0, col = "grey")
abline(h = ll.sens$sens$estimate[ll.sens$sens$alpha == 0] - 1000, lty = 2, col = "black")
polygon(x = c(sign(ll.sens$sens$alpha)*ll.sens$sens$rsqs, rev(sign(ll.sens$sens$alpha)*ll.sens$sens$rsqs)),
        y = c(ll.sens$sens$lower, rev(ll.sens$sens$upper)),
        col = rgb(0.5, 0.5, 0.5, alpha = 0.5),
        border = NA)
points(x = ll.sens$partial.r2, y = rep(0, length(ll.sens$partial.r2)), pch = 4)

##dev.off()
